After many hours of exploration, here is what I found:

Because smooth.spline algorithm chooses spar instead of lambda, it is only possible to (sort of) fix spar. However, lambda is a function of spar and another variable matrix. So fixing spar does not fix lambda necessarily. I have not found an easy way to extract the smoother matrix out of smooth.spline. However, for the purpose of computing variance, the algorithm provided in https://stat.ethz.ch/pipermail/r-help/2006-June/108471.html (fix spar instead of df) is a close estimate of the true smoother matrix. The variance computed from $SyS^T$, where $S$ is the estimated smoother matrix, is pretty close to the one computed from the correct smoother matrix.

Another R package "assist" has a function "ssr()" that also does smooth spline regression. It is not as powerful as smooth.spline. But the built-in function "hat.ssr()" gives the true smoother matrix of the model obtained from "ssr()".